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13.  ABSTRACT  (Maximum  200  words) 


The  object  of  this  research  is  to  study  the  effects  of  surface  topography,  near-surface  (sedimentary)  struc¬ 
ture  and  the  associated  small-scale  heterogeneities  on  regional  wave  propagation  which  is  critical  for  both 
discrimination  and  yield  estimation  in  monitoring  a  Comprehensive  Test  Ban  Treaty  and  the  current  Nuclear 
Non-Proliferation  Treaty.  We  aim  at  developing  a  hybrid  method  which  couples  the  recently  developed  fast 
screen  propagator  theory  and  methods  (Wu,  1994;  Wu  and  Xie,  1994;  Wu  and  Huang,  1995)  with  a  modified 
Boundary  Integral  Equation  (BIE)  method  to  treat  the  influences  of  both  volume  heterogeneities  and  irregular 
interfaces,  including  the  influence  of  surface  topography.  We  adopt  Chen’s  Global  Generalized 
Reflection/Transmission  Matrix  (GGRTM)  method  as  the  major  element  in  our  hybrid  method.  As  the  first  step, 
we  test  both  the  generalized  screen  one-way  wave  method  and  the  GGRTM  method,  and  develop  the  connection 
of  the  two  algorithms  for  the  two-dimensional  SH  case.  Chen’s  theory  has  been  modified  for  this  purpose  and 
the  connection  formulas  have  been  derived  and  numerically  tested.  The  excellent  agreement  between  seismo¬ 
grams  for  direct  propagation  and  propagation  using  the  connection  formulas  proves  the  correctness  of  the  theory 
and  the  feasibility  of  the  methodology. 
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1  Introduction 


The  study  of  path  effects  of  complex  structure  and  heterogeneities  on  the 
excitation  and  propagation  of  regional  phases  in  different  areas  remains  crit¬ 
ical  for  both  discrimination  and  yield  estimation  procedures  for  monitoring 
the  CTBT.  The  problem  will  be  more  severe  in  the  case  of  Non-Proliferation 
monitoring,  in  which  the  potential  nuclear  tests  may  occur  in  very  differ¬ 
ent  geological  and  geophysical  environments.  Today,  regional  waves  are  one 
of  the  most  important  indicators  for  monitoring  purpose.  Due  to  the  com¬ 
plexity  involved  in  the  regional  phase  propagation,  synthetic  simulation  will 
play  an  important  role  in  areas  that  lack  enough  data  to  rely  on  expirical 
methods.  To  meet  these  requirements,  the  ultimate  goal  is  to  develop  a 
computationally  viable  technique  for  calculating  high-frequency  (1-25  Hz) 
synthetic  seismograms  in  regional  distance  (  >  1000  km)  for  three  dimen¬ 
sional,  heterogeneous  (on  large  and  small  scales)  crustal  structures  including 
rough  surface  and  interfaces. 

In  the  past,  boundary  integral  equation  (BIE)  or  boundary  element  (BE) 
methods  have  been  extensively  used  to  study  the  effects  of  topography  or 
sedimentary  basin  structures  on  ground  motions  at  the  surface.  BE  has  been 
also  used  to  study  the  Lg  blockage  problem  with  limited  success.  Blockage  is 
assumed  to  be  caused  by  coastlines,  mountains  and  sudden  change  of  crustal 
thickness.  However,  two-dimensional  numerical  simulations  of  blockage  by 
large-scale  crustal  structures  have  not  succeeded  in  matching  the  observa¬ 
tions  (Campillo  et  al.,  1993;  Gibson  and  Campillo,  1994).  Most  simulations 
are  either  for  surface  topography  or  for  irregular  structure  beneath  a  flat 
surface  (sedimentary  layer)  due  to  the  restriction  of  computational  complex¬ 
ity.  However,  the  combination  of  both  surface-topography  and  sedimentary 
structure  may  have  more  dramatic  influence.  An  irregular  surface  and  low- 
velocity  layer  can  both  trap  part  of  the  Lg  energy  into  the  surface  layer 
and  scatter  the  Lg  wave  out  of  the  crustal  waveguide.  Existing  methods 
are  also  not  capable  of  simulating  the  combined  effects  of  both  large-scale 
structure  and  the  associated  small-scale  heterogeneities.  Irregular  topogra¬ 
phy  and  near-surface  structure  are  the  manifestation  of  past  and/ or  present 
tectonic  processes  which  often  produce  crustal  heterogeneities  at  different 
scales.  The  effects  of  the  small-scale  (wavelength-scale)  heterogeneities  must 
be  taken  into  consideration  in  modeling  blockage  and  other  Lg  propagation, 
scattering  and  attenuation  phenomena. 
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We  are  developing  a  new  hybrid  numerical  method  by  combining  the  gen¬ 
eralized  screen  method  with  the  boundary  integral  equation  method.  The 
generalized  screen  method  can  handle  wave  propagation  in  heterogeneous 
waveguide  with  modest  topography.  The  method  is  based  on  one-way  wave 
equation  theory  (Wu,  1994;  1996;  Wu  and  Xie,  1994;  Wu  and  Huang,  1995). 
In  the  crustal  waveguide  environment,  major  wave  energy  is  carried  by  for¬ 
ward  propagating  waves,  including  forward  scattered  waves,  and  therefore 
the  neglect  of  backscattered  waves  in  the  modeling  will  not  change  the  main 
features  of  regional  phases  in  most  cases.  By  neglecting  backscattering  in  the 
theory,  the  method  becomes  a  forward  marching  algorithm  in  which  the  next 
step  propagation  depends  only  on  the  present  value  of  wavefield  in  a  trans¬ 
verse  cross-section  and  the  heterogeneities  between  the  two  cross-sections. 
The  saving  of  computing  time  and  storage  is  enormous.  This  makes  it  a  very 
efficient  method  and  can  propagate  high  frequency  regional  signals  to  very 
long  distances. 

Modest  surface  topography  can  be  modeled  by  coordinate  transformation 
in  the  generalized  screen  method.  The  algorithm  for  handling  the  topography 
is  still  in  the  process  of  development.  On  the  other  hand,  the  boundary  inte¬ 
gral  equation  method  has  the  flexibility  to  incorporate  complex  topographic 
features  into  the  model.  However,  since  matrix  operations  are  involved,  the 
boundary  integral  equation  method  is  not  efficient.  When  the  ratio  of  model 
dimension  to  wavelength  is  too  large,  the  computation  time  and  memory  re¬ 
quirement  become  formidable.  This  problem  can  be  circumvented  through 
a  hybrid  method.  The  hybrid  method  will  combine  the  advantages  of  the 
above  mentioned  two  methods  and  avoid  their  disadvantages.  The  Lg  phases 
generated  by  the  source  are  propagated  to  a  certain  distance  with  the  gen¬ 
eralized  screen  method.  Then,  the  output  will  be  used  as  the  input  to  the 
boundary  integral  equation  method,  and  the  later  is  used  to  calculate  the  in¬ 
teraction  between  Lg  wave  and  the  complex  waveguide  structure  with  rough 
topographic  features.  This  approach  provide  us  a  possibility  to  investigate 
the  interaction  between  Lg  wave  and  crustal  waveguides  having  complicated 
structures  including  severe  topography  for  long  distance  propagation. 

As  the  first  year’s  effort,  the  screen  method  has  been  successfully  devel¬ 
oped  for  a  crustal  waveguide  for  2-D  SH-wave  propagation.  The  boundary 
integral  equation  method  has  been  tested  and  connected  to  a  laterally  ho¬ 
mogeneous  crust  model  for  2-D  SH-wave  propagation.  Numerical  examples 
showed  the  feasibility  of  this  approach.  The  next  year’s  work  will  be  to  put 
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the  topographic  features  into  the  hybrid  method,  and  test  the  method  for 
some  realistic  models,  such  as  those  for  paths  across  the  Tibet  Plateau. 


2  Generalized  Screen  Method 

For  an  isotropic  2D  elastic  medium,  the  SH  and  the  P-SV  waves  are  decou¬ 
pled.  Here  we  treat  only  the  SH  problem  to  demonstrate  the  applicability 
of  the  screen  propagators  to  crustal  waveguide.  Under  such  a  circumstance, 
the  equation  of  motion  becomes 

-<.V(r)u(r)  =  (1) 

where  oj  is  the  frequency,  r  =  (a;,  z)  is  a  2D  position  vector,  u  is  transverse 
displacement,  p  is  the  density  of  the  medium,  and  p,  is  the  shear  rigidity. 
We  decompose  the  parameters  of  the  elastic  medium  and  the  total  wave  field 
into 

P  =  Po  + 
p  =  Po  Sp 

u  =  u^  +  U  (2) 

where  po  and  po  are  parameters  of  the  background  medium,  8p  and  6p  are 
corresponding  perturbations,  u°  is  the  primary  field  and  U  is  the  scattered 
field.  Then,  the  SH  wave  equation  can  be  rewritten  as 

PqV^U  +  oj'^poU  =  —[u^6pu  +  V  •  SpS/u]  ,  (3) 

or 

(V^  +  e)U{r)  =  -eF{r)u{r) ,  (4) 

where  k  =  ojfv  is  the  wavenumber  in  the  background  medium  and  v  is  the 
background  S  wave  velocity  defined  by 

V  =  VPo/po  (5) 

In  the  right-hand  side  of  (4),  F{r)  is  a  perturbation  operator 

Fit)  =  e,(t)  +  pV .  ,  (6) 
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with 


/>0 

(7) 

(8) 

Equation  (4)  is  a  scalar  Helmholtz  equation.  With  a  half-space  scalar 
Green’s  function  the  scattered  field  U  can  be  written  as 


f7(ri)  =  j^dhg’^{ri,r)F{r)u{r)  , 


(9) 


where  the  2D  integration  is  over  the  volume  V  including  all  the  hetero¬ 
geneities  in  the  modeling  space.  Under  the  forward-scattering  approxima¬ 
tion,  the  total  field  and  Green’s  function  under  the  integration  in  above 
equation  can  be  replaced  by  their  forward-scattering  approximated  counter¬ 
parts,  and  the  field  can  be  calculated  by  a  one-way  marching  algorithm  along 
the  x-direction  using  a  dual  domain  technique. 


2.1  Wide-angle  Screen  Approximation 

The  half-space  model  can  be  sliced  into  thin-slabs  perpendicular  to  the  prop¬ 
agation  direction.  Weak  scattering  condition  holds  for  each  thin-slab.  For 
each  forward  step,  the  forward-scattered  field  by  the  thin-slab  is  calculated 
and  added  to  the  primary  field  so  that  the  updated  field  becomes  the  incident 
field  for  the  next  thin-slab.  The  formulas  of  the  dual-domain  implementation 
are  summarized  as  follows: 


U{xi,lQ  =  (10) 

where 

U,{xr,K,)  =  ik  r  dxe^'^^^^-^'>C[-Sp{z)uo{z)]  (11) 

Jx'  ^ 

U^{xi,K,)  =  ik  j  "  {c[e ^,{z)dxUo{z)]  -  i5[^£^(^)5^uo(^)]| 

^  (12) 
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where  C[f{z)]  and  <?[/(2)]  are  the  cosine  and  sine  transforms,  defined  by 

roo 

^[/(^)]  =  /  dz2cos{K;,z)f{z) 

Jo 

Jroo 

<  dz2sm{K,z)f{z)  (13) 

0 

In  Eq.  (11)  and  (12),  uq,  BxUq  and  BzUq  can  be  calculated  by 

ZTT  J— oo 

=  (14) 

and 

BxUoix,z)  = 

BzUoix,z)  =  iS-^[e^^'^^-^''>^uo{x\Ki)]  (15) 

Eq.  (11),  (12),  (14)  and  (15)  are  the  dual-domain  expressions  of  the  wide- 

angle  screen  propagator  for  half-space  SH  problems. 

The  procedure  can  be  summarized  as  follows. 

1.  Cosine  transform  the  incident  fields  at  the  entrance  of  each  thin-slab 
into  wavenumber  domain. 

2.  Free  propagate  in  wavenumber  domain  and  calculate  the  primary  field 
and  its  gradient  within  the  slab. 

3.  At  each  horizontal  position  within  the  slab,  inverse  cosine/sine  trans¬ 
form  the  primary  field  and  its  gradients  into  space  domain  and,  then  interact 
with  the  medium  perturbations  Sp  and  e^. 

4.  Cosine/sine  transform  the  distorted  fields  into  wavenumber  domain 
and  perform  the  divergence  operations  to  get  the  scattered  fields 

5.  Calculate  the  primary  field  at  the  slab  exit  and  add  to  the  scattered 
field  to  form  the  total  field  as  the  incident  field  at  the  entrance  of  the  next 
thin-slab. 

6.  Continue  the  procedure  iteratively. 
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2.2  Small  Angle  Screen  Approximation  and  the  Phase- 
Screen  Propagator 

When  the  energy  of  crustal  guided  waves  are  carried  mainly  by  small-angle 
waves  (with  respect  to  the  horizontal  direction),  the  small  angle  approxi¬ 
mations  can  be  invoked  to  simplify  the  theory  and  calculations.  Under  the 
phase-screen  approximation,  the  heterogeneous  half-space  is  represented  by  a 
series  of  half-screens  embedded  in  the  homogeneous  background  half-space. 
The  wave  propagates  between  screens  in  the  wavenumber  domain  and  in¬ 
teracts  with  phase-screens  in  the  space  domain.  The  interaction  is  only 
a  phase-delay  operator  (multiplication  in  space  domain).  The  formula  for 
dual-domain  implementation  is 

it(iri ,  uq(^xi j  Ji.z}  T  U (3^1,  ) 

/‘CO 

_  j  dz2cos{Kzz)[l  +  ik2Ss{z)]uo{x',z) 

Jo 

^  (16) 

where  exp[2ikSs{z)]  is  the  phase  delay  operator.  The  procedure  can  be  sum¬ 
marized  as  follows. 

1.  Cosine  transform  the  incident  field  at  the  starting  plane  into  wavenum¬ 
ber  domain  and  free  propagate  to  the  screen. 

2.  Inverse  cosine  transform  the  incident  field  into  space  domain  and 
interact  with  the  shear  slowness  screen  (phase-screen)  to  get  the  transmitted 
field. 

3.  Cosine  transform  the  transmitted  field  into  wavenumber  domain  and 
free  propagate  to  the  next  screen. 

4.  Repeat  the  propagation  and  interaction  screen-by-screen  to  the  bound¬ 
ary  of  the  model  space. 

2.3  Treatment  of  the  Moho  Discontinuity 

The  Moho  discontinuity  can  be  treated  in  two  ways.  One  is  to  put  the 
impedance  boundary  conditions  in  the  formulation,  the  other  is  to  treat  the 
parameter  changes  as  perturbations  and  therefore  be  incorporated  into  the 
screen  interaction.  The  former  has  the  advantage  of  computational  efficiency. 
The  latter  has  the  flexibility  of  handling  irregular  interfaces.  Here,  we  adopt 
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the  latter  approach  and  check  the  validity  of  the  perturbation  approach  for 
the  Moho  discontinuity  by  a  reflectivity  method  and  a  finite  difference  algo¬ 
rithm. 


3  Global  Generalized  Reflection  Transmis¬ 
sion  Matrix  Method 

The  discretization  of  BIE  can  be  done  by  integration  of  the  Green’s  func¬ 
tion  either  in  space  domain  (e.g.  Sanches-Sesma  and  Campillo,  1991),  or  in 
wavenumber  domain  using  the  discrete  wave  number  representation  (Bou- 
chon,  1985;  Campillo  and  Bouchon,  1985;  Chen,  1990,  1995,  1996).  In  the 
latter  approach,  the  singularity  problem  of  the  Green’s  function  is  avoided 
by  using  truncated  series.  The  wavenumber  domain  BIE  has  another  ad¬ 
vantage  that  it  can  be  easily  extended  to  the  case  of  multilayered  media 
with  irregular  interfaces.  In  Bouchon  et  al.  (1989),  propagator  matrices 
are  used  to  relate  equivalent  force  distributions  on  neighboring  interfaces. 
Chen  (1990,  1995,  1996)  related  the  fields  at  neighboring  layers  by  global 
reflection/transmission  coefficients  and  then  derived  the  global  generalized 
R/T  coefficients  to  relate  observations  and  sources.  In  these  methods,  the 
dimensionality  of  the  linear  system  to  solve  are  independent  of  the  num¬ 
ber  of  layers  involved.  The  computation  time  increases  only  linearly  with 
the  number  of  interfaces.  For  this  reason,  we  adopt  Chen’s  GGRTM  (Global 
Generalized  Reflection/Transmission  Matrix)  method  as  the  candidate  in  our 
hybrid  method. 

The  GGRTM  can  be  viewed  as  an  extension  of  the  reflectivity  method 
for  horizontally  layered  case  to  an  irregularly  layered  case,  and  it  has  been 
demonstrated  to  be  an  accurate  and  effective  method  to  simulate  seismic 
waves  in  laterally  varying  layered  media  (Chen,  1991,  1995,  1996).  For  ex¬ 
ample,  for  the  scattering  problem  due  to  a  semi-circular  canyon  (shown  in 
Figure  1),  GGRTM  can  provide  very  accurate  results.  Figures  2  and  3  show 
the  comparisons  of  the  results  (solid  lines)  computed  by  GGRTM  with  the 
analytical  solutions  of  Trifunac  (dotted  lines)  for  various  normalized  frequen¬ 
cies,  showing  excellent  agreement  between  them.  It  is  known  that  in  this 
semi-circular  canyon  model,  there  are  two  sharp  edges.  Many  other  meth¬ 
ods,  e.g.,  Aki-Larner  method,  T-matrix  method  and  other  high-frequency 
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asymptotic  methods,  fail  to  provide  correct  solutions. 

3.1  Connection  Formulation 

Assume  domain  II  is  the  model  space  we  are  interested  in  and  the  field  in 
domain  I  is  easy  to  calculate  by  other  less  expensive  methods.  According  to 
the  representation  theorem,  wave-fields  inside  domain  II  can  be  expressed  as 

+CO  ✓  ^ 

=  J  I  r^(x')  +  u^{x')n{z')— 

0  ^ 

Where  and  are  the  displacement  and  traction  fields  on  the  vertical 
boundary  surface  dividing  domain  I  and  II,  and  can  be  calculated  using 
methods  valid  in  domain  I,  //  is  the  shear  rigidity,  and  is  the  Green’s 
function  in  domain  II  which  will  be  calculated  by  GGRTM. 

3.2  Algorithm  of  Computing  Synthetic  Lg  Waves 

Having  the  connection  formulation,  we  can  use  GGRTM  to  compute  syn¬ 
thetic  Lg  waves.  The  step-by-step  procedure  of  applying  GGRTM  to  com¬ 
puting  a  synthetic  seismogram  in  a  general  irregularly  layered  medium  can 
be  summarized  as  follows. 

Step  1 

Calculate  the  interface  matrices  for  each  interface,  Qn  ’  ’ 

P|^|\  Pj^^l  ;  for  j=l,2,  ...,  N,  by  carrying  out  the  integrals  over  each 
interface.  These  interface  matrices  contain  the  structural  information  of  the 
media  and  are  defined  as  (Chen,  1990) 


"^G^\x,x!)dz'  (17) 


-1 


2ui^'>L 


Ll2 

j  +  4'^^}  exp[i2['!f^(a;,  n,  m)]dx  , 


-L/2 


(18) 


L/2 

j  -  i/(-’)}exp[iS[-’|^(a:,n,m)]dx 


-L/2 


(19) 
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Figure  1:  The  configuration  of  the  scattering  problem  due  to  a  semi-circular 
canyon  and  an  incident  plane  wave,  where  a  is  the  radius  of  the  canyon,  and 
0  is  the  angle  of  incident  wave. 


L/2 

(^n)„  =  2~^  -f  -  i/^^^}exp[zS}f(a;,n,m)]da;  ,  (20) 


exp[iS}f  (a;, n, m)]dx  ,  (21) 


-L/2 

L/2 


2v^^>L 


=  /  {l  +  ^^(^)P}^^exp[zS5f(a:,n,m)]da:  ,  (22) 


-L/2 

L/2 


2l/n^L 


/  {l  +  ''H^)r}^^exp[iSjf(a:,n,m)]dx  ,  (23) 

2l/n 

L/2 

y  {l  +  [<^^^)(x)]^y%xp[^HJ^^(x,n,?7^)]o^x  ,  (24) 


-L/2 


L/2 


/p0)\  _  ^  rn 

V^TWn  “ 


(j)\  — 


and 

(pS>). 


-L/2 

L/2 


2d^)u^J^L 


J  {l  +  [^^^^3:)]^}  exp[2S^//(x,7?.,m)]d3;  ,  (25) 


-L/2 


9 


0, 


-3.  -2.  -1.  0.  1. 

distance  (x/a) 


— j 

3. 


77=0,25 


Figure  2:  The  frequency  responses  of  a  semi-circular  canyon  to  vertical  inci¬ 
dent  SH-wave  for  various  normalized  frequencies.  The  solid  lines  denote  our 
results  and  the  dotted  lines  denote  the  exact  solutions. 
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2. 


J7  =  1.0  6=20° 


Figure  3:  The  same  response  as  Figure  2,  except  that  the  incident  angle  is 
30^ 

where  is  the  height  of  the  topography  for  the  jth  interface,  and 

kn  =  2TrnL  -  {k-aY  ,  and  >  0, 

Sjf  (a:,n,m)  =  {k^  -  kr,)x  +  , 


s['^/(x,n,m)  =  (fern  -  kn)x  +  ^^(x)  -  ^ 

and 

3[f  (x,n,m)  =  {km  -  kn)x  -  |A^(^"^)(x)|  ; 

for  j  =  1,  2,  ...,  N. 

Where 
Step  2 

Calculate  the  global  modified  reflection  and/or  transmission  matrices,  { 
},  from  the  interface  matrices  using  the  following  formulas: 


■  13  (i)  rp(i)  " 

np(i)  T>U) 

L-i-u  ^Ti  -I 


Wtt 

_Q(i+i) 

Hii 


T>i^) 

^TT 

p(i+0 

'*^11 


iQiy" 


p(-?) 

p(i+i) 


-1 


E(^I 


0+1) 


■^mm 


(26) 


U 


and 


R-U  =-Qu  KtT’ESa,  (27) 

Where  e2„  and  EW,  are  diagonal  matrices  given  by 

E^n  =  diagonal  {exp[zn(-'^(^2n  “  ±1.  ±2, ...}  , 

and 

E^L  =  diagonal  n  =  0,  ±1,  ±2, ...} . 

These  global  modified  reflection  and/or  transmission  matrices  describe 
the  reflection  and/or  transmission  effects  due  to  single  interface  regardless  of 
the  influences  from  other  existing  interfaces. 

Step  3 

Compute  the  global  generalized  reflection  and/or  transmission  matrices, 
T|p  and  from  the  global  modified  reflection  and/or  trans¬ 
mission  matrices  through  following  recursive  formulas: 

p  (°)  _  p  (‘2) 

t;?  =  [I-RSfRi',-‘>l-'TS?,  /or;  =  l,2,...,iV;  (28) 

’dU)  _  rpOO-p  (i”l)rp(i) 

and 

=  0 

Tif  =  /ori  =  iV,iV-l,...,2,l.  (29) 

p(j)  _  p(i)  I  rp(j)-pO'+l)r|-i(j) 

These  global  generalized  reflection  and/or  transmission  matrices  represent 
the  total  reflection  and/or  transmissions  due  to  the  multi-irregular  layers. 
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Step  4 

Compute  the  expansion  coefficients  of  displacement  spectra  at  free  sur¬ 
face, 

by  using  the  formula 


^(o)  /a(1)  ,  T(1)-(2)  ,  rp(l)r^(2)^(3)  ,  , 

^  j  ^nun\S|  +J-TT®T  +  ■‘■TT  ■'■TT®T  +  -  +  ... X^l 


(l)rp(2)  ^(iV).(N+l)J  ^ 

(30) 


where  is  the  equivalent  source  term  for  the  jth  layer  derived  by  the 
representation  theorem  and 


2l/n^L 


i</'  =  {l  -  Rir'’Rl?}"  (sf  +  .  (31) 

LI2  40)(x) 

J  dx  j  +  -^^J]dz  , 


-L/2  ^0-i)(x) 
L/2  «(»(x) 


(32) 


”  2i/«  L_ii^ 

(33) 

for  j  =1,  2,  ...,  N,  N+1; 

and  ,  .  .  /)  1 

f^^\x^z)  =  <.T^^\x,z)  +  .  (34) 

Step  5 

Calculate  the  displacement  spectra  at  the  free  surface  by  using  the  fol¬ 
lowing  formula; 


M 

VF('’^[a:,^(°^(x),a;]  =  ^  exp  |A^(°^(x)|}  .  (35) 

^=-M 


Taking  the  Fourier  transform  on  the  above  frequency  domain  solution,  we 
can  finally  obtain  the  time  domain  solution,  i.e.,  the  synthetic  seismogram. 
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4  Numerical  Simulations 

4.1  Numerical  Simulations  for  Screen  Method 

In  this  section,  we  give  examples  of  using  the  half-space  phase-screen  algo¬ 
rithm  for  regional  wave  propagation.  First,  in  Figure  4  we  show  the  accuracy 
of  the  method  by  comparing  the  synthetic  seismograms  generated  by  the 
screen  method  (thick  lines)  with  those  calculated  by  a  reflectivity  method 
(thin  lines)  for  a  flat  crustal  model.  The  crust  has  a  thickness  of  32  km 
and  a  shear  wave  velocity  of  3.5  km/s.  The  mantle  beneath  the  crust  has 
a  shear  velocity  of  4.5  km/s.  The  source  function  is  a  Ricker  wavelet  with 
a  dominant  frequency  of  1.0  Hz.  Except  for  near  vertical  reflections,  where 
one  way  wave  equation  methods  have  difficulty,  the  results  show  excellent 
agreement.  For  long  distance  regional  waves,  the  contribution  of  near  ver¬ 
tical  reflections  is  negligible.  Next,  we  show  the  accuracy  of  the  method 
by  comparing  synthetic  seismograms  generated  by  this  method  with  those 
generated  by  a  finite  difference  algorithm  (Xie  and  Lay,  1994).  For  the  finite- 
difference  method,  a  fourth-order  elastic  SH-wave  code  is  used  to  calculate 
the  synthetic  seismograms.  The  spatial  sampling  interval  is  0.125  km  and  the 
time  interval  is  0.015  second.  For  the  screen  method,  the  spatial  sampling 
interval  is  0.25  km  in  vertical  direction  and  the  screen  interval  is  1.0  km. 
A  Gaussian  derivative  is  used  as  the  source  time  function  for  both  methods. 
Because  of  the  computational  intensity  of  the  finite  difference  method,  we  did 
the  comparison  at  short  propagation  distances.  Shown  on  the  top  of  Figure 

5  is  the  crust  model  used  to  calculate  synthetic  seismograms;  on  the  bottom, 
synthetic  seismograms  along  a  vertical  profile  at  an  epicenter  distance  of  250 
km.  The  thin  lines  are  from  the  finite  difference  method  and  the  thick  lines 
are  from  the  generalized  screen  method.  The  source  is  located  at  a  depth  of 
2  km.  Excellent  agreement  can  be  seen. 

Figure  6  shows  the  snap  shots  from  the  Screen  method  at  50  sec.  for  flat, 
narrowing  and  broadening  crustal  waveguides  (from  top  to  bottom,  respec¬ 
tively).  The  source  is  located  at  the  top-left  corner  at  depth  2  km.  The  de¬ 
velopment  of  mantle  wave  and  head  wave,  and  the  formation  of  crust  guided 
wave  as  multiple  reflections  between  the  free  surface  and  Moho  discontinuity 
can  be  clearly  seen.  For  the  inhomogeneous  models,  wave  diffraction,  leakage 
to  the  mantle,  wavefront  distortion  and  increase  of  wavefleld  complexity  can 
be  also  seen  clearly.  From  the  comparison  it  is  seen  that  the  passage  of  narrow 
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SP(thick  line)  compared  with  Reflectivity(thin  line) 


Figure  4:  Comparison  of  synthetic  seismograms  along  the  surface  calculated 
by  the  screen  method  and  reflectivity  method  for  a  flat  crustal  model  (32  km 
thick).  The  source  function  is  a  Ricker  wavelet  with  dominant  frequency  of 
LO  Hz. 


Vertical  Distance  (km) 


SP(thick  solid  line)  &  FD(thin  line)  for  inhomogeneous  model 


Figure  5:  Comparison  of  synthetic  seismograms  along  a  vertical  profile  at 
the  distance  of  250  km  calculated  by  the  screen  method  (thick  lines)  and 
a  finite-difference  method  (thin  lines)  for  a  laterally  varying  crustal  model 
shown  on  the  top  panel. 
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Coiiparison  of  Uaue  Propagation 
in  Uarious  Crustal  Uaue  Guides 


(  t  =  50  sec) 


Figure  6:  Snap  shots  at  50  sec.  for  various  crustal  waveguides. 
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crustal  waveguide  has  greater  effect  on  Lg  leakage  than  the  broad  passage. 
In  the  latter  case  although  the  wavefronts  are  complicated  due  to  scattering 
at  the  edges,  most  of  the  energy  is  still  trapped  in  the  crust,  different  from 
the  case  of  a  narrow  passage  in  which  a  large  percentage  of  energy  leaks  into 
the  mantle.  This  example  demonstrates  the  potential  of  the  method  as  a 
tool  for  investigating  the  path  effects  of  different  crustal  structures. 

Figure  7  shows  the  synthetic  seismograms  by  the  screen  method  for  the 
Flora- Asnes  crust  model  in  the  NORSAR  region.  The  parameters  for  this 
model  are  listed  in  Table  1.  The  model  has  a  low  velocity  top  layer  1  km 
thick  and  a  velocity  discontinuity  at  depth  15  km.  The  receivers  are  on  the 
surface.  A  Ricker  wavelet  is  used  for  this  simulation  with  /o  =  1  Hz.  Shown 
in  the  upper  panel  are  short  distances  (up  to  350  km);  the  lower  panel  shows 
long  distances  (up  to  1000  km).  In  this  case  the  Lg  group  is  formed  by 
multiple  reflections  of  the  Moho  and  the  crustal  discontinuities,  complicated 
by  the  low  velocity  sedimental  layer. 


Table  1:  Flora- Asnes  crust  model 


thickness  {km) 

Vs  [kmls) 

P  (S'/cm^) 

1.00 

3.00 

2.60 

3.46 

2.80 

3.76 

3.00 

infinity 

4.65 

3.30 

The  following  example  shows  the  potential  capability  of  this  method  for 
long  distance  high-frequency  synthetic  seismograms  in  a  laterally  varying 
structure.  Figure  8  shows  the  laterally  varying  crust  model  used  in  the 
calculation.  Figure  9  shows  the  high-frequency  synthetic  seismograms  on 
the  surface  at  distances  up  to  1000  km.  The  center  frequency  is  5  Hz  with 
the  maximum  frequency  of  10  Hz.  In  comparison,  the  low-frequency  (/c  = 
1  Hz,  fmax  =  2  Hz)  synthetic  seismograms  are  shown  in  Figure  10.  It  is 
clear  that  without  high-frequency  content,  many  of  the  distinctive  features 
associated  with  Lg  measurements  can  not  be  adequately  modeled.  In  other 
words,  a  proper  simulation  method  with  the  capability  to  generate  accurate 
high-frequency  signals  is  a  necessity  for  the  purpose  of  investigating  regional 
phases.  The  generalized  screen  method  with  its  high  efficiency  serves  well  as 
an  important  element  of  the  hybrid  method. 


Figure  7:  synthetic  seismograms  for  the  Flora- Asnes  crust  model  in  the  NORSAR 
region.  The  parameters  for  this  model  are  listed  in  Table  1. 
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►eismogram  of  Flora-Asnes 


Parameters  of  Crustal  Model 


Layer 

Vsfkm/sec) 

DensityCg/cm ) 

Thickness(km) 

1 

3.00 

2.60 

1.00 

2 

3.46 

2.80 

14.00 

3 

3.76 

3.00 

22.00 

4 

4.65 

3.30 

Half-Space 

Crustal  Model 

Figure  8:  An  inhomogeneous  crustal  model  used  in  the  calculation  of  h-f 
synthetic  seismograms.  Shown  in  the  upper  panel  are  model  parameters  and 
the  lower  panel  gives  the  geometry  of  the  model.  The  receivers  are  on  the 
surface  and  shown  by  triangles. 
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Synthetic  seismograms  (fc=5.0Hz,fmaxs5l  O.OHz),  Time(Sec) 


Distance(km) 


Distance(km) 


Figure  10;  low-frequency  (/c  =  1  Hz,  fmax  =  2  Hz)  synthetic  seismograms 
on  the  surface  at  distances  up  to  1000  km  for  an  inhomogeneous  crustal 
waveguide  (Figure  8). 


4.2  Numerical  Test  for  BIE  Method 

To  test  the  validity  of  our  hybrid  method,  we  consider  a  trivial  case:  a 
laterally  homogeneous  layered  model.  This  problem  can  be  fully  solved  by 
reflectivity  method.  To  test  our  algorithm,  we  use  our  hybrid  method  to  syn¬ 
thesize  the  seismograms,  then  check  the  results  with  the  reflectivity  method. 
The  test  model  is  a  single  layer  crustal  model.  The  velocities  and  densi¬ 
ties  of  the  crust  and  mantle  are  3.5  fcm/sec,  2.8  g/cm^,  4.5  kmfsec  and  3.2 
glcm^  ,  respectively.  The  thickness  of  the  crust  is  32  km,  seismic  source 
is  buried  at  2^=2  km  and  X5=0  km.  Receiver  is  placed  at  zo=0  km  and 
a;o=250  km.  The  connection  boundary  is  located  at  a:=150  km.  The  syn¬ 
thetic  seismogram  of  reflectivity  method  is  plotted  in  Figure  11a  and  the 
synthetic  seismogram  from  GGRTM  is  shown  in  Figure  11b.  Comparison  of 
these  two  seismograms  shows  an  excellent  agreement,  confirming  the  validity 
of  the  connection  scheme  for  our  hybrid  method. 

The  computer  code  for  calculating  general  irregular  media  is  under  de¬ 
velopment  at  this  stage,  and  expected  to  be  finished  soon.  We  will  then 
calculate  synthetic  Lg  waves  propagating  through  an  arbitrarily  irregular 
layered  medium  to  study  the  influence  of  surface  topography  and  interface 
irregularities. 

5  Conclusion  and  Discussions 

We  have  derived  the  connection  formulas  for  our  hybrid  scheme  and  its  va¬ 
lidity  has  been  proved  by  numerical  tests.  Both  generalized  screen  method 
and  boundary  integral  equation  method  have  been  tested  for  the  waveguide 
environment.  The  algorithm  for  seismogram  synthesis  in  arbitrarily  irregular 
layered  media  is  under  development.  We  will  also  study  the  approximation 
involved  in  the  so-called  Rayleigh  Ansatz  method,  or  the  Aki-Larner  method 
(Aki  and  Lamer,  1970).  Aki-Larner  method  is  a  wavenumber  domain  im¬ 
plemented  and  approximated  BIE  method.  It  is  much  faster  than  the  strict 
BIE  method  and  therefore  can  simulate  large  3D  topography  and  interface 
problems.  Horike  et  al.  (1990)  has  applied  the  method  to  non-axisymmetric 
3D  surface  structure  problems.  We  will  test  the  accuracy  and  speed  of  AL 
method  by  comparing  it  with  the  strict  BIE  method  (such  as  Chen’s  GGRTM 
method)  and  incorporate  the  approximation  into  our  hybrid  method. 
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displacement  displacement 


Comparison  of  Synthetic  Seismograms 


Figure  11:  Comparison  of  synthetic  seismograms  for  a  laterally  homogeneous 
layered  crustal  model.  A:  synthetic  seismogram  from  a  reflectivity  method, 
and  B:  synthetic  seismogram  from  the  hybrid  method. 
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